We create GLMMs to investigate whether the raising type or colony has an influence on a) the number of fledglings per potential mother (incl. breeding probability) and b) the number of fledglings per nest (only actual breeding events and incl. temporarily added females as mothers).

Setup

## for non-CRAN packages please keep install instruction
## but commented so it is not run each time, e.g.
# devtools::install_github("EcoDynIZW/template")

## libraries used in this script
## please add ALL LIBRARIES NEEDED HERE
## please remove libraries from the list that are not needed anymore 
## at a later stage
library("DHARMa")
## Registered S3 methods overwritten by 'lme4':
##   method                          from
##   cooks.distance.influence.merMod car 
##   influence.merMod                car 
##   dfbeta.influence.merMod         car 
##   dfbetas.influence.merMod        car
## This is DHARMa 0.3.3.0. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa') Note: Syntax of plotResiduals has changed in 0.3.0, see ?plotResiduals for details
library("effects")
## Loading required package: carData
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
library("ggeffects")
library("ggplot2")
library("here")
library("lme4")
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
library("pwr")
library("tidyverse")

Data

f_fledg_mother <- readRDS(file = here("output", "data-proc", "f_fledg_mother.rds")) # from script "05a_Reproduction_analysis"          # chicks_year 
f_fledg_mother <- droplevels(f_fledg_mother)
nest_numbers_fm <- readRDS(file = here("output", "data-proc", "nest_numbers_fm.rds")) # from script "04_Survival_analysis"        # NBI_breeding_all
nest_numbers_fm <- droplevels(nest_numbers_fm)

Female fledglings per potential mother per year

Based on the values used for RR Baseline

names(f_fledg_mother)
## [1] "parent_f"      "hatch_year"    "nr_f_fledg"    "raising_type"  "breeding_area"

GLMM of the number of fledglings per potential mother per raising type

Data exploration

# Plot the number of fledglings
ggplot(data = f_fledg_mother, mapping = aes(x = nr_f_fledg)) + 
  geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# Plot of the number of fledglings per raising type per year
f_fledg_mother %>%
  group_by(raising_type, hatch_year) %>%
  mutate(sum.fledg = sum(nr_f_fledg)) %>%
  ggplot(mapping = aes(x = hatch_year, y = sum.fledg)) + 
  geom_line(aes(group = raising_type, color = raising_type))

# Sum of the breeding events per raising type
table(f_fledg_mother$raising_type)
## 
## foster parent        parent 
##            47            16
# Mean per group
RT_year <- aggregate(x = f_fledg_mother, by = list(Raising = f_fledg_mother$raising_type, Year = f_fledg_mother$hatch_year), 
                     FUN = mean)
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
RT_mean <- aggregate(x = RT_year, by = list(RT_year$Raising), 
            FUN = mean)
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA

GLM

null_model_f_fledg_raising <- glm(formula = nr_f_fledg ~ 1, 
                                  data = f_fledg_mother, 
                                  family = "poisson")
summary(null_model_f_fledg_raising)
## 
## Call:
## glm(formula = nr_f_fledg ~ 1, family = "poisson", data = f_fledg_mother)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.0541  -1.0541  -1.0541   0.5354   2.2868  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.5878     0.1690  -3.478 0.000506 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 70.964  on 62  degrees of freedom
## Residual deviance: 70.964  on 62  degrees of freedom
## AIC: 128.63
## 
## Number of Fisher Scoring iterations: 5
f_fledg_raising_glm <- glm(formula = nr_f_fledg ~ raising_type,
                           data = f_fledg_mother,
                           family = "poisson")
summary(f_fledg_raising_glm) # non significant
## 
## Call:
## glm(formula = nr_f_fledg ~ raising_type, family = "poisson", 
##     data = f_fledg_mother)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.1299  -1.1299  -0.7906   0.4177   2.8628  
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)  
## (Intercept)         -0.4490     0.1826  -2.459   0.0139 *
## raising_typeparent  -0.7142     0.4830  -1.479   0.1393  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 70.964  on 62  degrees of freedom
## Residual deviance: 68.387  on 61  degrees of freedom
## AIC: 128.05
## 
## Number of Fisher Scoring iterations: 6

Plot

#dredge(raising_mod)
plot(allEffects(f_fledg_raising_glm),type='response')

GLMM

f_fledg_raising_glmm <- glmer(formula = nr_f_fledg ~ raising_type + (1|parent_f), 
                              data = f_fledg_mother, 
                              family = "poisson")
summary(f_fledg_raising_glmm) # not significant
## Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
##  Family: poisson  ( log )
## Formula: nr_f_fledg ~ raising_type + (1 | parent_f)
##    Data: f_fledg_mother
## 
##      AIC      BIC   logLik deviance df.resid 
##    125.3    131.7    -59.7    119.3       60 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.8567 -0.5852 -0.3750  0.3852  2.7081 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  parent_f (Intercept) 0.6684   0.8176  
## Number of obs: 63, groups:  parent_f, 31
## 
## Fixed effects:
##                    Estimate Std. Error z value Pr(>|z|)  
## (Intercept)         -0.8428     0.3509  -2.402   0.0163 *
## raising_typeparent  -0.9307     0.7060  -1.318   0.1874  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## rsng_typprn -0.142

Plot

f_fledg_raising_predi <- ggpredict(f_fledg_raising_glmm, terms = "raising_type")
f_fledg_raising_predi
## # Predicted counts of nr_f_fledg
## # x = raising_type
## 
## x             | Predicted |       95% CI
## ----------------------------------------
## foster parent |      0.43 | [0.22, 0.86]
## parent        |      0.17 | [0.04, 0.73]
## 
## Adjusted for:
## * parent_f = 0 (population-level)
plot(x = f_fledg_raising_predi, 
     show.title = F, 
     limits = c(0,1)) + 
  labs(x = "Raising type",
       y = "Number of female fledglings per potential mother", 
       tags = "") +
  theme_classic() + 
  theme( 
    axis.line = element_line(colour = "black"),
    axis.text = element_text(size=14),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank(), 
    plot.margin = unit(c(7.5, 12.5, 5.5, 15.5), "points"),
    axis.title.x = element_text(colour="black", size = 15),
    axis.title.y = element_text(colour="black",size=15, vjust = 5))

Testing the models

We use the DHARMa package.

Prepare it

f_fledg_raising_simu <- simulateResiduals(fittedModel = f_fledg_raising_glmm)
plot(f_fledg_raising_simu)

Dispersion

testDispersion(f_fledg_raising_simu)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.94162, p-value = 0.88
## alternative hypothesis: two.sided

Zero-inflation

testZeroInflation(f_fledg_raising_simu)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.9154, p-value = 0.488
## alternative hypothesis: two.sided

Temporal Autocorrelation

testTemporalAutocorrelation(f_fledg_raising_simu)
## DHARMa::testTemporalAutocorrelation - no time argument provided, using random times for each data point

## 
##  Durbin-Watson test
## 
## data:  simulationOutput$scaledResiduals ~ 1
## DW = 2.1776, p-value = 0.4782
## alternative hypothesis: true autocorrelation is not 0

GLMM of the number of fledglings per potential mother per colony

Data exploration

# Plot the number of fledglings
ggplot(data = f_fledg_mother, mapping = aes(x = nr_f_fledg)) + 
  geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# Plot of the number of fledglings per colony per year
f_fledg_mother %>%
  group_by(breeding_area, hatch_year) %>%
  mutate(sum.fledg = sum(nr_f_fledg)) %>%
  ggplot(mapping = aes(x = hatch_year, y = sum.fledg)) + 
  geom_line(aes(group = breeding_area, color = breeding_area))

# Sum of the breeding events per colony
table(f_fledg_mother$breeding_area)
## 
## Burghausen      Kuchl allocation 
##         28         28          7

Remove individuals with colony “allocation”

f_fledg_mother <- f_fledg_mother[f_fledg_mother$breeding_area != "allocation", ]
f_fledg_mother <- droplevels(f_fledg_mother)
table(f_fledg_mother$breeding_area) # half half
## 
## Burghausen      Kuchl 
##         28         28

GLM

null_model_f_fledg_col <- glm(formula = nr_f_fledg ~ 1, 
                              data = f_fledg_mother, 
                              family = "poisson")
summary(null_model_f_fledg_col)
## 
## Call:
## glm(formula = nr_f_fledg ~ 1, family = "poisson", data = f_fledg_mother)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.1180  -1.1180  -1.1180   0.4359   2.1591  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)   -0.470      0.169  -2.781  0.00542 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 62.719  on 55  degrees of freedom
## Residual deviance: 62.719  on 55  degrees of freedom
## AIC: 120.39
## 
## Number of Fisher Scoring iterations: 5
f_fledg_colony_glm <- glm(formula = nr_f_fledg ~ breeding_area,
                          data = f_fledg_mother,
                          family = "poisson")
summary(f_fledg_colony_glm) # not significant
## 
## Call:
## glm(formula = nr_f_fledg ~ breeding_area, family = "poisson", 
##     data = f_fledg_mother)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.1339  -1.1019  -1.1019   0.4116   2.1909  
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)  
## (Intercept)        -0.49899    0.24254  -2.057   0.0396 *
## breeding_areaKuchl  0.05716    0.33820   0.169   0.8658  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 62.719  on 55  degrees of freedom
## Residual deviance: 62.691  on 54  degrees of freedom
## AIC: 122.36
## 
## Number of Fisher Scoring iterations: 6

Plot

#dredge(raising_mod)
plot(allEffects(f_fledg_colony_glm),type='response')

GLMM

f_fledg_colony_glmm <- glmer(formula = nr_f_fledg ~ breeding_area + (1|parent_f), 
                             data = f_fledg_mother,  
                             family = "poisson")
summary(f_fledg_colony_glmm)
## Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
##  Family: poisson  ( log )
## Formula: nr_f_fledg ~ breeding_area + (1 | parent_f)
##    Data: f_fledg_mother
## 
##      AIC      BIC   logLik deviance df.resid 
##    121.1    127.2    -57.6    115.1       53 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.9485 -0.5909 -0.5461  0.5472  2.2144 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  parent_f (Intercept) 0.4809   0.6934  
## Number of obs: 56, groups:  parent_f, 25
## 
## Fixed effects:
##                    Estimate Std. Error z value Pr(>|z|)  
## (Intercept)        -0.77960    0.38967  -2.001   0.0454 *
## breeding_areaKuchl  0.02262    0.47364   0.048   0.9619  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## brdng_rKchl -0.606

Plot

f_fledg_colony_predi <- ggpredict(f_fledg_colony_glmm, terms = "breeding_area")
f_fledg_colony_predi
## # Predicted counts of nr_f_fledg
## # x = breeding_area
## 
## x          | Predicted |       95% CI
## -------------------------------------
## Burghausen |      0.46 | [0.21, 0.98]
## Kuchl      |      0.47 | [0.22, 1.01]
## 
## Adjusted for:
## * parent_f = 0 (population-level)
plot(x = f_fledg_colony_predi, 
     show.title = F, 
     limits = c(0,1.15)) + 
  labs(x = "Colony",
       y = "Number of female fledglings per potential mother", 
       tags = "") +
  theme_classic() + 
  theme( 
    axis.line = element_line(colour = "black"),
    axis.text = element_text(size=14),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank(), 
    plot.margin = unit(c(7.5, 12.5, 5.5, 15.5), "points"),
    axis.title.x = element_text(colour="black", size = 15),
    axis.title.y = element_text(colour="black",size=15, vjust = 5))

Testing the models

We use the DHARMa package.

Prepare it

f_fledg_colony_simu <- simulateResiduals(fittedModel = f_fledg_colony_glmm)
plot(f_fledg_colony_simu)

Dispersion

testDispersion(f_fledg_colony_simu)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.96286, p-value = 0.968
## alternative hypothesis: two.sided

Zero-inflation

testZeroInflation(f_fledg_colony_simu)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.92394, p-value = 0.608
## alternative hypothesis: two.sided

Temporal Autocorrelation

testTemporalAutocorrelation(f_fledg_colony_simu)
## DHARMa::testTemporalAutocorrelation - no time argument provided, using random times for each data point

## 
##  Durbin-Watson test
## 
## data:  simulationOutput$scaledResiduals ~ 1
## DW = 1.5061, p-value = 0.05923
## alternative hypothesis: true autocorrelation is not 0

Female and male fledglings per nest per year

Based on the values for RRNest ## GLMM of the number of fledglings per nest per raising type

names(nest_numbers_fm)
##  [1] "b_pair_date"       "b_parent_f"        "b_parent_f_id"     "b_parent_m"        "b_parent_m_id"     "b_breeding_colony" "b_count_eggs"      "b_count_hatching" 
##  [9] "b_count_fledgling" "breed_year"        "raising_type"

Data exploration

# Plot the number of fledglings
ggplot(data = nest_numbers_fm, mapping = aes(x = b_count_fledgling)) + 
  geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# Plot of the number of fledglings per raising type per year
nest_numbers_fm %>% 
  group_by(raising_type, breed_year) %>%
  mutate(sum.fledg = sum(b_count_fledgling)) %>%
  ggplot(mapping = aes(x = breed_year, y = sum.fledg)) + 
  geom_line(aes(group = raising_type, color = raising_type))

# Sum of the breeding events with reproduction per nest and per raising type
table(nest_numbers_fm$raising_type)
## 
## foster parent        parent         added 
##            31             3            32

Analysis with temporarily added females as mothers

GLM

null_model_fm_nest_raising_full <- glm(formula = b_count_fledgling ~ 1, 
                                       data = nest_numbers_fm, 
                                       family = "poisson")
summary(null_model_fm_nest_raising_full)
## 
## Call:
## glm(formula = b_count_fledgling ~ 1, family = "poisson", data = nest_numbers_fm)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.06706  -0.86863  -0.09432   0.55657   1.66615  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.75911    0.08422   9.014   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 72.857  on 65  degrees of freedom
## Residual deviance: 72.857  on 65  degrees of freedom
## AIC: 227.22
## 
## Number of Fisher Scoring iterations: 5
fm_nest_raising_full_glm <- glm(formula = b_count_fledgling ~ raising_type,
                                data = nest_numbers_fm,
                                family = "poisson")
summary(fm_nest_raising_full_glm) # non significant
## 
## Call:
## glm(formula = b_count_fledgling ~ raising_type, family = "poisson", 
##     data = nest_numbers_fm)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.30940  -0.82416  -0.04514   0.60937   1.27256  
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         0.72490    0.12500   5.799 6.66e-09 ***
## raising_typeparent  0.25593    0.37500   0.682    0.495    
## raising_typeadded   0.04347    0.17354   0.251    0.802    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 72.857  on 65  degrees of freedom
## Residual deviance: 72.409  on 63  degrees of freedom
## AIC: 230.77
## 
## Number of Fisher Scoring iterations: 5

Plot

#dredge(raising_mod)
plot(allEffects(fm_nest_raising_full_glm), type='response')

GLMM

fm_nest_raising_full_glmm <- glmer(formula = b_count_fledgling ~ raising_type + (1|b_parent_f), 
                                   data = nest_numbers_fm, 
                                   family = "poisson")
## boundary (singular) fit: see ?isSingular
summary(fm_nest_raising_full_glmm) # not significant
## Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
##  Family: poisson  ( log )
## Formula: b_count_fledgling ~ raising_type + (1 | b_parent_f)
##    Data: nest_numbers_fm
## 
##      AIC      BIC   logLik deviance df.resid 
##    232.8    241.5   -112.4    224.8       62 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6330 -0.7409 -0.0449  0.6511  1.4289 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  b_parent_f (Intercept) 0        0       
## Number of obs: 66, groups:  b_parent_f, 21
## 
## Fixed effects:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         0.72490    0.12500   5.799 6.66e-09 ***
## raising_typeparent  0.25593    0.37500   0.682    0.495    
## raising_typeadded   0.04347    0.17354   0.251    0.802    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) rsng_typp
## rsng_typprn -0.333          
## rsng_typddd -0.720  0.240   
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular

Plot

fm_nest_raising_full_predi <- ggpredict(fm_nest_raising_full_glmm, terms = "raising_type")
fm_nest_raising_full_predi
## # Predicted counts of b_count_fledgling
## # x = raising_type
## 
## x             | Predicted |       95% CI
## ----------------------------------------
## foster parent |      2.06 | [1.62, 2.64]
## parent        |      2.67 | [1.33, 5.33]
## added         |      2.16 | [1.70, 2.73]
## 
## Adjusted for:
## * b_parent_f = 0 (population-level)
plot(x = fm_nest_raising_full_predi, 
     show.title = F, 
     limits = c(0,6)) + 
  labs(x = "Raising type",
       y = "Number of fledglings per nest (m + f)", 
       tags = "") +
  theme_classic() + 
  theme( 
    axis.line = element_line(colour = "black"),
    axis.text = element_text(size=14),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank(), 
    plot.margin = unit(c(7.5, 12.5, 5.5, 15.5), "points"),
    axis.title.x = element_text(colour="black", size = 15),
    axis.title.y = element_text(colour="black",size=15, vjust = 5))

Testing the models

We use the DHARMa package.

Prepare it

fm_nest_raising_full_simu <- simulateResiduals(fittedModel = fm_nest_raising_full_glmm)
plot(fm_nest_raising_full_simu)

Dispersion

testDispersion(fm_nest_raising_full_simu)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.89881, p-value = 0.272
## alternative hypothesis: two.sided

Zero-inflation

testZeroInflation(fm_nest_raising_full_simu)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.375, p-value = 0.32
## alternative hypothesis: two.sided

Temporal Autocorrelation

testTemporalAutocorrelation(fm_nest_raising_full_simu)
## DHARMa::testTemporalAutocorrelation - no time argument provided, using random times for each data point

## 
##  Durbin-Watson test
## 
## data:  simulationOutput$scaledResiduals ~ 1
## DW = 1.9305, p-value = 0.7767
## alternative hypothesis: true autocorrelation is not 0

Analysis without temporarily added females as mothers

Here we consider only the raising types foster parent and biological parent

nest_numbers_fm_bpfp <- nest_numbers_fm[nest_numbers_fm$raising_type != "added",]

GLM

null_model_fm_nest_raising <- glm(formula = b_count_fledgling ~ 1, 
                                  data = nest_numbers_fm_bpfp, 
                                  family = "poisson")
summary(null_model_fm_nest_raising)
## 
## Call:
## glm(formula = b_count_fledgling ~ 1, family = "poisson", data = nest_numbers_fm_bpfp)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.05798  -0.66326  -0.08161   0.57021   1.68125  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.7503     0.1179   6.367 1.93e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 32.558  on 33  degrees of freedom
## Residual deviance: 32.558  on 33  degrees of freedom
## AIC: 114.34
## 
## Number of Fisher Scoring iterations: 5
fm_nest_raising_glm <- glm(formula = b_count_fledgling ~ raising_type,
                           data = nest_numbers_fm_bpfp,
                           family = "poisson")
summary(fm_nest_raising_glm) # not significant
## 
## Call:
## glm(formula = b_count_fledgling ~ raising_type, family = "poisson", 
##     data = nest_numbers_fm_bpfp)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.30940  -0.62940  -0.04514   0.60937   1.27256  
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          0.7249     0.1250   5.799 6.66e-09 ***
## raising_typeparent   0.2559     0.3750   0.682    0.495    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 32.558  on 33  degrees of freedom
## Residual deviance: 32.122  on 32  degrees of freedom
## AIC: 115.9
## 
## Number of Fisher Scoring iterations: 5

Plot

#dredge(raising_mod)
plot(allEffects(fm_nest_raising_glm), type='response')

GLMM

fm_nest_raising_glmm <- glmer(formula = b_count_fledgling ~ raising_type + (1|b_parent_f), 
                                   data = nest_numbers_fm_bpfp, 
                                   family = "poisson")
## boundary (singular) fit: see ?isSingular
summary(fm_nest_raising_glmm) # not significant
## Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
##  Family: poisson  ( log )
## Formula: b_count_fledgling ~ raising_type + (1 | b_parent_f)
##    Data: nest_numbers_fm_bpfp
## 
##      AIC      BIC   logLik deviance df.resid 
##    117.9    122.5    -55.9    111.9       31 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6330 -0.5669 -0.0449  0.6511  1.4289 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  b_parent_f (Intercept) 0        0       
## Number of obs: 34, groups:  b_parent_f, 14
## 
## Fixed effects:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          0.7249     0.1250   5.799 6.66e-09 ***
## raising_typeparent   0.2559     0.3750   0.682    0.495    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## rsng_typprn -0.333
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular

Plot

fm_nest_raising_predi <- ggpredict(fm_nest_raising_glmm, terms = "raising_type")
fm_nest_raising_predi
## # Predicted counts of b_count_fledgling
## # x = raising_type
## 
## x             | Predicted |       95% CI
## ----------------------------------------
## foster parent |      2.06 | [1.62, 2.64]
## parent        |      2.67 | [1.33, 5.33]
## 
## Adjusted for:
## * b_parent_f = 0 (population-level)
plot(x = fm_nest_raising_predi, 
     show.title = F, 
     limits = c(0,6)) + 
  labs(x = "Raising type",
       y = "Number of fledglings per nest (m + f)", 
       tags = "") +
  theme_classic() + 
  theme( 
    axis.line = element_line(colour = "black"),
    axis.text = element_text(size=14),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank(), 
    plot.margin = unit(c(7.5, 12.5, 5.5, 15.5), "points"),
    axis.title.x = element_text(colour="black", size = 15),
    axis.title.y = element_text(colour="black",size=15, vjust = 5))

Testing the models

Testing with the DHARMa package.

Prepare it

fm_nest_raising_simu <- simulateResiduals(fittedModel = fm_nest_raising_glmm)
plot(fm_nest_raising_full_simu)

Dispersion

testDispersion(fm_nest_raising_simu)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.85891, p-value = 0.312
## alternative hypothesis: two.sided

Zero-inflation

testZeroInflation(fm_nest_raising_simu)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.2054, p-value = 0.816
## alternative hypothesis: two.sided

Temporal Autocorrelation

testTemporalAutocorrelation(fm_nest_raising_simu)
## DHARMa::testTemporalAutocorrelation - no time argument provided, using random times for each data point

## 
##  Durbin-Watson test
## 
## data:  simulationOutput$scaledResiduals ~ 1
## DW = 1.5941, p-value = 0.2277
## alternative hypothesis: true autocorrelation is not 0

GLMM of the number of fledglings per nest per colony

Data exploration

# Plot the number of fledglings
ggplot(data = nest_numbers_fm, mapping = aes(x = b_count_fledgling)) + 
  geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# Plot of the number of fledglings per raising type per year
nest_numbers_fm %>% 
  group_by(b_breeding_colony, breed_year) %>%
  mutate(sum.fledg = sum(b_count_fledgling)) %>%
  ggplot(mapping = aes(x = breed_year, y = sum.fledg)) + 
  geom_line(aes(group = b_breeding_colony, color = b_breeding_colony))

# Sum of the breeding events with reproduction per nest and per raising type
table(nest_numbers_fm$b_breeding_colony)
## 
## Burghausen      Kuchl 
##         40         26

GLM

null_model_fm_nest_colony <- glm(formula = b_count_fledgling ~ 1, 
                                 data = nest_numbers_fm, 
                                 family = "poisson")
summary(null_model_fm_nest_colony)
## 
## Call:
## glm(formula = b_count_fledgling ~ 1, family = "poisson", data = nest_numbers_fm)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.06706  -0.86863  -0.09432   0.55657   1.66615  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.75911    0.08422   9.014   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 72.857  on 65  degrees of freedom
## Residual deviance: 72.857  on 65  degrees of freedom
## AIC: 227.22
## 
## Number of Fisher Scoring iterations: 5
fm_nest_colony_glm <- glm(formula = b_count_fledgling ~ b_breeding_colony,
                          data = nest_numbers_fm,
                          family = "poisson")
summary(fm_nest_colony_glm) # non significant
## 
## Call:
## glm(formula = b_count_fledgling ~ b_breeding_colony, family = "poisson", 
##     data = nest_numbers_fm)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.13037  -0.81506  -0.03506   0.62019   1.73666  
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              0.7178     0.1104   6.500 8.02e-11 ***
## b_breeding_colonyKuchl   0.1016     0.1707   0.595    0.552    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 72.857  on 65  degrees of freedom
## Residual deviance: 72.504  on 64  degrees of freedom
## AIC: 228.86
## 
## Number of Fisher Scoring iterations: 5

Plot

#dredge(raising_mod)
plot(allEffects(fm_nest_colony_glm), type='response')

GLMM

fm_nest_colony_glmm <- glmer(formula = b_count_fledgling ~ b_breeding_colony + (1|b_parent_f), 
                             data = nest_numbers_fm, 
                             family = "poisson")
## boundary (singular) fit: see ?isSingular
summary(fm_nest_colony_glmm) # not significant
## Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
##  Family: poisson  ( log )
## Formula: b_count_fledgling ~ b_breeding_colony + (1 | b_parent_f)
##    Data: nest_numbers_fm
## 
##      AIC      BIC   logLik deviance df.resid 
##    230.9    237.4   -112.4    224.9       63 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.50640 -0.73335 -0.03492  0.66351  2.06037 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  b_parent_f (Intercept) 0        0       
## Number of obs: 66, groups:  b_parent_f, 21
## 
## Fixed effects:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              0.7178     0.1104   6.500 8.02e-11 ***
## b_breeding_colonyKuchl   0.1016     0.1707   0.595    0.552    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## b_brdng_clK -0.647
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular

Plot

fm_nest_colony_predi <- ggpredict(fm_nest_colony_glmm, terms = "b_breeding_colony")
fm_nest_colony_predi
## # Predicted counts of b_count_fledgling
## # x = b_breeding_colony
## 
## x          | Predicted |       95% CI
## -------------------------------------
## Burghausen |      2.05 | [1.65, 2.55]
## Kuchl      |      2.27 | [1.76, 2.93]
## 
## Adjusted for:
## * b_parent_f = 0 (population-level)
plot(x = fm_nest_colony_predi, 
     show.title = F, 
     limits = c(0,3)) + 
  labs(x = "Breeding colony",
       y = "Number of fledglings per nest (m + f)", 
       tags = "") +
  theme_classic() + 
  theme( 
    axis.line = element_line(colour = "black"),
    axis.text = element_text(size=14),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank(), 
    plot.margin = unit(c(7.5, 12.5, 5.5, 15.5), "points"),
    axis.title.x = element_text(colour="black", size = 15),
    axis.title.y = element_text(colour="black",size=15, vjust = 5))

Testing the models

Testing with the DHARMa package.

Prepare it

fm_nest_colony_simu <- simulateResiduals(fittedModel = fm_nest_colony_glmm)
plot(fm_nest_colony_simu)

Dispersion

testDispersion(fm_nest_colony_simu)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.89981, p-value = 0.304
## alternative hypothesis: two.sided

Zero-inflation

testZeroInflation(fm_nest_colony_simu)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.3798, p-value = 0.32
## alternative hypothesis: two.sided

Temporal Autocorrelation

testTemporalAutocorrelation(fm_nest_colony_simu)
## DHARMa::testTemporalAutocorrelation - no time argument provided, using random times for each data point

## 
##  Durbin-Watson test
## 
## data:  simulationOutput$scaledResiduals ~ 1
## DW = 2.1516, p-value = 0.5357
## alternative hypothesis: true autocorrelation is not 0

Session Info

## DO NOT REMOVE!
## We store the settings of your computer and the current versions of the
## packages used to allow for reproducibility
Sys.time()
## [1] "2022-01-06 20:04:48 CET"
#git2r::repository() ## uncomment if you are using GitHub
sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22000)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252    LC_MONETARY=German_Germany.1252 LC_NUMERIC=C                    LC_TIME=German_Germany.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] pwr_1.3-0         lme4_1.1-25       Matrix_1.2-18     ggeffects_1.0.0   effects_4.2-0     carData_3.0-4     DHARMa_0.3.3.0    lubridate_1.7.9   timelineS_0.1.1  
## [10] viridisLite_0.3.0 forcats_0.5.0     stringr_1.4.0     dplyr_1.0.2       purrr_0.3.4       readr_1.4.0       tidyr_1.1.2       tibble_3.0.3      tidyverse_1.3.0  
## [19] survminer_0.4.8   ggpubr_0.4.0      ggplot2_3.3.2     survival_3.2-7    here_0.1         
## 
## loaded via a namespace (and not attached):
##   [1] minqa_1.2.4       colorspace_1.4-1  ggsignif_0.6.0    ellipsis_0.3.1    rio_0.5.16        sjlabelled_1.1.7  showtext_0.9      rprojroot_1.3-2   snakecase_0.11.0 
##  [10] estimability_1.3  fs_1.5.0          rstudioapi_0.11   farver_2.0.3      showtextdb_3.0    remotes_2.2.0     fansi_0.4.1       xml2_1.3.2        codetools_0.2-16 
##  [19] splines_4.0.3     knitr_1.30        pkgload_1.1.0     jsonlite_1.7.1    nloptr_1.2.2.2    broom_0.7.4       km.ci_0.5-2       dbplyr_1.4.4      compiler_4.0.3   
##  [28] d6_0.1.0.0        httr_1.4.2        backports_1.1.10  assertthat_0.2.1  survey_4.0        cli_2.0.2         htmltools_0.5.0   prettyunits_1.1.1 tools_4.0.3      
##  [37] gtable_0.3.0      glue_1.4.2        tinytex_0.26      Rcpp_1.0.5        cellranger_1.1.0  vctrs_0.3.4       nlme_3.1-149      iterators_1.0.13  lmtest_0.9-38    
##  [46] insight_0.11.1    xfun_0.18         ps_1.4.0          openxlsx_4.2.3    testthat_2.3.2    rvest_0.3.6       lifecycle_0.2.0   devtools_2.3.2    statmod_1.4.35   
##  [55] rstatix_0.6.0     MASS_7.3-53       zoo_1.8-8         scales_1.1.1      hms_0.5.3         yaml_2.2.1        curl_4.3          memoise_1.1.0     gridExtra_2.3    
##  [64] KMsurv_0.1-5      stringi_1.5.3     desc_1.2.0        gap_1.2.2         foreach_1.5.1     boot_1.3-25       pkgbuild_1.1.0    zip_2.1.1         rlang_0.4.7      
##  [73] pkgconfig_2.0.3   evaluate_0.14     lattice_0.20-41   labeling_0.3      processx_3.4.4    tidyselect_1.1.0  magrittr_1.5      bookdown_0.20     R6_2.4.1         
##  [82] generics_0.0.2    DBI_1.1.0         pillar_1.4.6      haven_2.3.1       foreign_0.8-80    withr_2.3.0       nnet_7.3-14       abind_1.4-5       modelr_0.1.8     
##  [91] crayon_1.3.4      car_3.0-10        survMisc_0.5.5    rmarkdown_2.4     sysfonts_0.8.1    usethis_1.6.3     grid_4.0.3        readxl_1.3.1      data.table_1.13.2
## [100] blob_1.2.1        callr_3.5.0       rmdformats_0.3.7  reprex_0.3.0      digest_0.6.25     xtable_1.8-4      munsell_0.5.0     mitools_2.4       sessioninfo_1.1.1